function Y = featuresPseudoZernike(Img, featureArgs,ve)
%FEATURESPSEUDOZERNIKE spocita pseudozernikovy momenty
%                 Img - vstupni fotka
%                 ve - elipsa kolem obliceje [r1 r2 cx cy theta]
%                 featureArgs - vektor obsahujici polozky n1 a n2 -
%                                n1 - minimalni stupen
%                                n2 - maximalni stupen
%                 Y - vektor (n2+n1+1)(n2-n1+2)/2 zernikovych momentu


Img = im2double(Img);
[SizeY, SizeX] = size(Img);

%nastaveni defaultnich hodnot parametru
if nargin < 2
    n1 = 2;
    n2 = 7;
else
    n1 = featureArgs(1);
    n2 = featureArgs(2);
end;

if nargin < 3
   [x0, y0] = findCenter(Img);
   r = 0.8*min([(SizeY - y0) y0 (SizeX - x0) x0]);
else
   x0 = ve(3);
   y0 = ve(4);
   r = min(ve([1 2]));
end  


%vypocet vzdalenosti (D) a uhlu (A) 
[X Y] = meshgrid((1:SizeX) - x0, (1:SizeY) - y0);
D = sqrt(X .^2 + Y .^2) ./ r;
A = atan(Y ./ X);

%prevod dat na vektory 
D = reshape(D, SizeX*SizeY, 1);
A = reshape(A, SizeX*SizeY, 1);
Img = reshape(Img, SizeX*SizeY, 1);

%vynechani bodu mimo kruh
In = D < 1;
D = D(In);
A = A(In);
Img = Img(In);
cs = length(D);

%vypocet sinu a kosinu
SinA = zeros(cs, 1+n2);
CosA = zeros(cs, 1+n2);
CosA(:, 1) = 1;
for i = 1:n2
    SinA(:,1+i) = sin(i*A);
    CosA(:,1+i) = cos(i*A);
end

%vynasobeni obrazem
SinA = SinA .* repmat(Img, 1, n2+1);
CosA = CosA .* repmat(Img, 1, n2+1);

%vypocet pseudo-zernikovych polynomu
R = cell(n2+1,n2+1);
R{1+0, 1+0} = ones(cs, 1);
for q = 1:n2
    R{1+q, 1+q} = D .* R{q, q};
    R{1+q, q} = (2*q+1)*R{1+q,1+q} + (2*q)*R{q,q};
end
if n2 >= 2
    D2 = R{3,3};
    for q = 0:(n2-2)
        for p = (q+2):n2
            L1 = ((2*p+1)*(2*p))/((p+q+1)*(p-q));
            L2 = -2*p+((p+q)*(p-q-1)*L1)/(2*p-1);
            L3 = (2*p-1)*(p-1)-0.5*(p+q-1)*(p-q-2)*L1+2*(p-1)*L2;
            R{1+p,1+q} = (L1*D2+L2).*R{p,1+q}+L3*R{p-1,1+q};
        end
    end
end

%vypocet momentu
ys = (n2+n1+1)*(n2-n1+2)/2;
Y = zeros(ys, 1);
i = 1;
for p = n1:n2
    for q = 0:p
        Y(i) = (p+1)/cs * sum(R{1+p,1+q} .* SinA(:, 1+q));
        Y(i+1) = (p+1)/cs * sum(R{1+p,1+q} .* CosA(:, 1+q));
        i = i+2;
    end
end

end

